load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl" ; Dev = "png" ; myDate = 20160516 ;NOMADS = "http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global"+myDate+"/rtofs_glo_2ds_forecast_daily_prog" NOMADS = "http://nomads.ncep.noaa.gov:9090/dods/rtofs/rtofs_global"+myDate+"/rtofs_glo_2ds_forecast_3hrly_prog" filename = NOMADS ; = url print("Fetching URL: ") print(" "+filename) ;#// + "/rtofs.20160125/rtofs_glo_2ds_n001_1hrly_prog.nc" exists = isfilepresent(filename) if(.not.exists) then print("OPeNDAP isfilepresent test unsuccessful.") print("Either file doesn't exist, or NCL does not have OPeNDAP capabilities on this system") exit else print("OPeNDAP isfilepresent test successful.") f = addfile(filename,"r") vars = getfilevarnames(f) print(vars) ; should be (in any order): if(.not.any(ismissing(vars))) then do i=0,dimsizes(vars)-1 printFileVarSummary (f,vars(i)) end do end if end if ; Define boundary of data ; Assign lat/lon boundaries to the data ; Longitudes between 0 and 74 E need to have 360 added to them ; myNlat = 32.0 mySlat = 10.0 myWlon = 260.0 myElon = 300.0 ;myNlat = 26.0 ;mySlat = 21.0 ;myWlon = 260.0 ;myElon = 300.0 ;myNlat = 90.0 ;mySlat = -90.0 ;myWlon = 74.0 ;myElon = 434.0 ; Get lat/lon/time info first lat = f->lat print ("Done reading lat") lon = f->lon print ("Done reading lon") time = f->time(0) print ("Done reading time") ; ; Use just one time when "0" mytime = 0 ;mytime = ":" Nlat = closest_val(myNlat, lat) print ("Nlat = " + Nlat + " Lat = " + lat(Nlat) ) Slat = closest_val(mySlat, lat) print ("Slat = " + Slat + " Lat = " + lat(Slat) ) Wlon = closest_val(myWlon, lon) print("Wlon = " + Wlon + " Lon = " + lon(Wlon) ) Elon = closest_val(myElon, lon) print("Elon = " + Elon + " Lon = " + lon(Elon) ) ; ;exit ; mylat = lat(Slat:Nlat) mylon = lon(Wlon:Elon) printVarSummary (mylat) printVarSummary (mylon) exit delete(lat) delete(lon) lat = f->lat(Slat:Nlat) lon = f->lon(Wlon:Elon) print ("Done reading subset lat/lon") printVarSummary(lat) printVarSummary(lon) sst = f->sst(mytime,:,Slat:Nlat,Wlon:Elon) ;sst = f->sst print ("Done reading sst") printVarSummary(sst) print (sst) u_velocity = f->u_velocity(mytime,:,Slat:Nlat,Wlon:Elon) ;u_velocity = f->u_velocity print ("Done reading U") printVarSummary(u_velocity) print (u_velocity) v_velocity = f->v_velocity(mytime,:,Slat:Nlat,Wlon:Elon) ;v_velocity = f->v_velocity print ("Done reading V") printVarSummary(v_velocity) print (v_velocity) speed = sqrt((u_velocity*u_velocity) + (v_velocity*v_velocity)) print ("Done calculating Speed") print (speed) exit ;;;;;; ;;;;;; ; Finished fetching data ; loop over all times ; This will do all times in the file ;do i=0,dimsizes(time)-1 do i=0,0 print ("time(i) = "+time(i)) ; Convert to UTC time. ; Array to hold month abbreviations. Don't store anything in index ; '0' (i.e. let index 1=Jan, 2=Feb, ..., index 12=Dec). ; month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \ "Oct","Nov","Dec"/) ; utc_date = cd_calendar(time(i), 0) ; ; Store return information into more meaningful variables. ; year = tointeger(utc_date(:,0)) ; Convert to integer for month = tointeger(utc_date(:,1)) ; use sprinti day = tointeger(utc_date(:,2)) hour = tointeger(utc_date(:,3)) minute = tointeger(utc_date(:,4)) second = utc_date(:,5) ; ; Write out strings in the format "hhZ dd mmm yyyy". ; date_str = sprinti("%0.2iZ ", hour) + sprinti("%0.2i ", day) + \ month_abbr(month) + " " + sprinti("%0.4i", year) fileDateFormat = "%Y%N%D_%H%M" fileDate = cd_string(time(i),fileDateFormat) ; print(" ===== Working on "+date_str+" =====") plotVarName = (/ "SST", "Speed" /) do plotLoop = 0, 1 ; Device Type is set at the beginning of the script DevTitle = "Images2/" + "cuba_" + plotVarName(plotLoop)+ "_" + fileDate ;***************************************************** ; Create plot ;***************************************************** wks = gsn_open_wks(Dev,DevTitle) ; select output ; gsn_define_colormap(wks,"BlWhRe") ; choose colormap gsn_define_colormap(wks,"rainbow") ; choose colormap icolor = NhlNewColor(wks,1.0,.89,.76) ; add bisque to colormap res = True ; plot mods desired ; Regional settings here res@gsnAddCyclic = False ; regional data res@mpGridSpacingF = 1.0 ; Lat/Lon grid spacing ; res@mpGridLatSpacingF = 0.1 ; res@mpGridLonSpacingF = 0.1 res@mpGridAndLimbOn = True res@mpGridLineThicknessF = 0.5 ; res@cnFillOn = True ; turn on color res@cnLinesOn = False res@gsnSpreadColors = True ; use full color map res@gsnSpreadColorEnd = -3 ; don't use added bisque res@lbLabelAutoStride = True ; nice label bar labels res@lbOrientation = "Vertical" res@mpDataBaseVersion = "MediumRes" ; "LowRes" "HighRes" res@mpLandFillColor = "bisque" ; choose land color res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries res@mpUSStateLineColor = "red" ; make them red res@mpProjection = "CylindricalEquidistant" ;;;;;;; was here ; zoom in on region res@mpLimitMode = "LatLon" res@mpMinLatF = mySlat ; min(lat) res@mpMaxLatF = myNlat ; max(lat) res@mpMinLonF = myWlon ;min(lon) res@mpMaxLonF = myElon ;max(lon) res@mpCenterLonF = (myWlon+myElon)/2.0 res@mpCenterLatF = (mySlat+myNlat)/2.0 ; res@tmXMajorGrid = True ; res@tmYMajorGrid = True ; res@tmLabelAutoStride = True ; res@tmXBAutoPrecision = True ; res@tmXBBorderOn = True ; res@tmXBLabelsOn = True ; res@tmYLLabelsOn = True res@tmXBOn = True res@tmYLOn = True res@tmXBTickSpacingF = 1.0 res@tmYLTickSpacingF = 1.0 res@tmXBMinorOn = True ; res@tmXBMinorOutwardLengthF = 0.1 ; res@tmXBMinorPerMajor = 4 res@tiMainString = "Valid: " + date_str ; This helps to get the markers where we want them, on top. res@gsnFrame = False res@gsnDraw = False res@gsnMaximize = True ; Section for vectors res@vcRefMagnitudeF = 1.0 ; define vector ref mag res@vcRefLengthF = 0.045 ; define length of vec ref res@vcGlyphStyle = "CurlyVector" ; turn on curly vectors res@vcMinDistanceF = 0.017 ; thin vectors res@vcRefAnnoSide = "Right" ; Put ref vector under colorbar ; res@vcRefAnnoJust = "TopRight" ; Put ref vector under colorbar res@vcRefAnnoParallelPosF = 1.15 res@vcRefAnnoOrthogonalPosF = .05 ; move ref vector res@gsnScalarContour = True res@cnFillOn = True res@cnLinesOn = False ; turn off contour lines ; plot = gsn_csm_vector_scalar_map(wks,myU,myV,myFsst,res) ;;;;;;;;;;;;;;; ; For SST if (plotLoop .eq. 0) then ; res@cnLevelSelectionMode = "ExplicitLevels" ;"ManualLevels" res@cnLevelSelectionMode = "AutomaticLevels" ; res@cnLevelSelectionMode = "ManualLevels" ; manual contour levels res@cnLevels := fspan(73,82,37) ; res@cnMinLevelValF = 32.0 ; set min contour level ; res@cnMaxLevelValF = 96.0 ; set max contour level ; res@cnLevelSpacingF = 0.25 ; contour interval plot = gsn_csm_vector_scalar_map(wks,u_velocity(i,0,:,:),v_velocity(i,0,:,:),sst(i,0,:,:),res) end if if (plotLoop .eq. 1) then ; Do speed contours ;CGH; res@cnLevelSelectionMode = "ExplicitLevels" ;"ManualLevels" res@cnLevelSelectionMode = "AutomaticLevels" ; res@cnLevelSelectionMode = "ManualLevels" ; ; manual contour levels res@cnLevels := fspan(0.0,4.0,17) ; res@cnMinLevelValF = 0.0 ; set min contour level ;; res@cnMaxLevelValF = 96.0 ; set max contour level ; res@cnLevelSpacingF = 0.25 ; contour interval plot = gsn_csm_vector_scalar_map(wks,u_velocity(i,0,:,:),v_velocity(i,0,:,:),speed(i,0,:,:),res) end if ;;;;;;;;; ; plot = gsn_csm_contour_map(wks,myFsst,res) draw(plot) ; ; ;Charleston Finish at 32 37.077'N, 79 35.497'W ; 32.61795 280.40838 ; ; Ponce Inlet Start at 29 04.72'N, 80 53.70'W ; 29.07867 279.105 ; ; Marker gres = True gres@gsMarkerIndex = 15 gres@gsMarkerSizeF = 0.03 gres@gsMarkerThicknessF = 2.5 ;gres@gsMarkerColor = "white" ;gres@gsnFrame = False ;gres@gsnDraw = False ; marklat = (/ 32.61795, 29.07867 /) ; marklon = (/280.40838, 279.10500 /) marklat = (/ 23.2, 24.54 /) marklon = (/278.9, 278.19 /) marks = gsn_add_polymarker(wks,plot,marklon,marklat,gres) draw(plot) ; ; plotspd = gsn_csm_vector_scalar_map(wks,myU,myV,mySpeed,res) ; draw(plotspd) ; ; overlay(plot,plotspd) ; frame(wks) end do ; plotLoop end do ; all times